function [h1,l2,linf] = bberror(dofs,V,T,d,c,u,ux,uy)

[qw,qp] = quad_rule(20);

matd = vdm23(d, qp(:,1), qp(:,2), 1 - qp(:,1) - qp(:,2));

Id = diag(ones((d+1)*(d+2)/2,1));
desc = desc_pattern(d);
Du = d*de_cast_step(Id,d,1,0,-1,desc); % the direction derivate on v3-v1
Dv = d*de_cast_step(Id,d,0,1,-1,desc);  % the direction derivate on v3-v2
matu = vdm23(d-1, qp(:,1), qp(:,2), 1 - qp(:,1) - qp(:,2))*Du;
matv = vdm23(d-1, qp(:,1), qp(:,2), 1 - qp(:,1) - qp(:,2))*Dv;

nt = size(T,1); 
L2 = zeros(nt,1); H1 = zeros(nt,1); Linf = zeros(nt,1);
for k = 1:nt
    V1=V(T(k,1),:);V2=V(T(k,2),:);V3=V(T(k,3),:);
    x13 = V1(1)-V3(1);
    y13 = V1(2)-V3(2);
    x23 = V2(1)-V3(1);
    y23 = V2(2)-V3(2);
    J = x13*y23 - x23*y13; 
    
    qp_x = V3(1) + x13*qp(:,1) + x23*qp(:,2);
    qp_y = V3(2) + y13*qp(:,1) + y23*qp(:,2);

    fvec = feval(u, qp_x, qp_y);
    fxvec = feval(ux, qp_x, qp_y);
    fyvec = feval(uy, qp_x, qp_y);

    loc_c = c(dofs(:,k));
    fbb = matd*loc_c;
    fbbx= (matu*loc_c*y23 - matv*loc_c*y13)/J;
    fbby= (matv*loc_c*x13 - matu*loc_c*x23)/J;
    
    L2(k)=((fbb-fvec).^2)'*qw*J;
    Linf(k) = max(abs(fbb - fvec));
    H1(k)= ((fbby-fyvec).^2)'*qw*J+((fbbx-fxvec).^2)'*qw*J;
end
h1 = sqrt(sum(L2) + sum(H1));
l2 = sqrt(sum(L2));
linf = max(Linf);

return

function pattern = desc_pattern(d)
% this function may be execute only once in whole comput session.
% much like the function asce_pattern.
m = (d+1)*(d+2)/2;
pattern = zeros(m,3);
begin = 1;
for j = 0:d
    idx = (begin:(begin+j))';
    pattern(idx,1) = idx;
    pattern(idx,2) = idx + j + 1;
    pattern(idx,3) = idx + j + 2;
    begin = begin + j + 1;
end

return